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Abstract. We investigate the behavior of energy fluctuations in several models of granular gases maintained 
in a non-equilibrium steady state. In the case of a gas heated from a boundary, the inhomogeneities of 
C"j the system play a predominant role. Interpreting the total kinetic energy as a sum of independent but not 

identically distributed random variables, it is possible to compute the probability density function (pdf) 
"-^J , of the total energy. Neglecting correlations and using the analytical expression for the inhomogeneous 

. temperature profile obtained from the granular hydrodynamic equations, we recover results which have 

been previously observed numerically and which had been attributed to the presence of correlations. In 
O , order to separate the effects of spatial inhomogeneities from those ascribable to velocity correlations, we 

have also considered two models of homogeneously thermostated gases: in this framework it is possible to 
reveal the presence of non-trivial effects due to velocity correlations between particles. Such correlations 
stem from the inelasticity of collisions. Moreover, the observation that the pdf of the total energy tends to 
ly-^ ■ a Gaussian in the large system limit, suggests that they are also due to the finite size of the system. 

cn : 

, PACS. 45.70.-n Granular systems - 05.40.-a Fluctuation phenomena, random processes, noise, and Brow- 

nian motion - 47.57.Gc Granular flow 

in : 

O 1 Introduction 

! 

^3 In non-equilibrium statistical mechanics the emergence of non-Gaussian distributions is a remarkable feature that 
C5 ' marks an essential difference with typical results of equilibrium statistical mechanics. In particular, non-Gaussian 
L | behaviors for global quantities (i.e. system averaged) have been observed in different fields of physics p^, unveiling 
unexpected analogies between turbulent flows, equilibrium critical phenomena 0, and non-equilibrium instabilities 
[3] . Averaged quantities are expected to be free of the microscopic details of the system under study, in such a way that 
, their probability distribution function (pdf) should only depend on few parameters. Furthermore, they will contain 
relevant informations about the correlations among different parts of the system. 

Granular gases represent a simple example of many-particles system out of equilibrium. The simplest model consists 
in an assembly of smooth hard particles undergoing binary inelastic collisions. This model, known as the Inelastic Hard 
5-h Spheres (IHS) model, has stimulated a great interest for more than ten years. Its close similarity with the elastic hard 
9^ i spheres fluid has made possible the exploitation of the long-past known standard methods of kinetic theory. However 
the inclusion of energy dissipation in the collision law is sufficient to generate many surprising features that have been 
observed both in experiments and simulations, such as clustering, convection, non-equipartition of energy, etc. |4I5I6| . 

In this paper we study the pdfs of the total kinetic energy of granular gases maintained in a non-equilibrium steady 
state by an external driving mechanism which compensates for the energy loss due to collisions. We first focus on 
boundary driven gases, since this kind of driving mechanism can easily be achieved experimentally, vibrating at high 
frequencies the container of the gas or one of its walls |7I8I9I10| . Energy is hence provided at the boundaries and 
dissipated by inelastic collisions in the bulk of the system, leading to strong spatial inhomogeneities. A coarse-grained 
description of this steady state, in terms of density and temperature profiles, can be obtained by means of hydrodynamic 
equations |llj , which are a priori valid when the mean-free-path is small compared to the typical length scale of the 
gradients. The presence of non-uniform density and temperature profiles implies that the total energy is a sum of 
variables which are not identically distributed and, in principle, correlated. We will first investigate the consequences of 
the lack of homogeneity by assuming total independence among the velocities, i.e. absence of correlations. Interestingly 
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enough, this approximate approach leads to good estimates for the shape of the pdfs. Similarly, Bcrtin in a recent 
paper |12| showed that considering the sum of independent, but not identically distributed, random variables can be 
a useful way to explain (at least on a qualitative level) the emergence of non-Gaussian pdfs in other systems. 

We will also discuss another important characteristic of such inhomogeneous systems, i.e. the fact that they are 
often non-extensive. Since the energy is injected from the boundaries and the dissipation happens inside the volume, 
it is clear that increasing the volume at constant density will not a priori increase the total energy linearly with the 
total number of particles N. The hydrodynamic approach allows to show that the relevant variable is the number of 
layers of particles at rest ^Jj • 

In order to disentangle the effects of inhomogeneities from those of correlations, we focus, in the second part of 
the paper, on the energy fluctuations in homogeneous granular gases. We will enforce homogeneity by considering the 
system in a pure velocity space, i.e. eliminating any spatial degree of freedom. When inhomogeneities disappear, the 
effects that may contribute, in the thermodynamic limit, to a non-Gaussian pdf of the total energy are the correlations 
among velocities. Nonetheless such a behavior generally shows up when the system is near some critical point, i.e. when 
some correlation length (in velocity space) diverges. A general theory for correlation functions for the hard sphere fluid 
has been first derived in an d fruitfully applied to the homogeneous cooling state (HCS) for a gas of inelastic hard 
spheres [14j . This last result was obtained exploiting the fact that the dynamics is dominated by the hydrodynamic 
(i.e. slow) modes, which are exactly computable for the particular case of the HCS. In our homogeneous models the 
total (kinetic) energy is the sum of N identically distributed random variables. We will show, by numerical simulations, 
that in this state the central limit theorem applies and the energy pdf is a Gaussian in the large N limit, which means 
that the correlation length in the velocity space is finite. The interesting issue then concerns the computation of the 
variance of the rescaled pdf. This quantity depends on the one-particle and (with a pre- factor 1/N) two-particles 
distribution functions. It is deeply related to the finite size correlations between the velocities of the particles. In both 
cases of a stochastic and of a deterministic thermostat, its measure reveals the presence of non-trivial correlations. 
Remarkably, it only depends upon the restitution coefficient, but not upon other parameters such as the number of 
particles or the amplitude of the thermostating force. In the deterministic thermostated case, it is moreover in almost 
perfect agreement with the theoretical computations given in |14j for a gas of inelastic hard spheres. 

The structure of the paper is the following: in section II we discuss the gas driven at the boundaries, while in 
section III we consider the homogeneous thermostats. Conclusions are drawn in section IV. 



2 The inhomogeneous boundary driven gas 

The focus is here on the energy fluctuations of a granular gas in the case where energy is injected through a bound- 
ary, typically a vibrating wall. This kind of system corresponds to a realistic experimental setup |7l8l9ll()j and has 
therefore been widely studied both numerically and analytically [1511611711811 llliJ| . As an important consequence of 
the boundary heating, the density and the temperature are not homogeneous: a heat flux is present, which does not 
verify the Fourier law. This feature is well described by kinetic theory and in good agreement with the hydrodynamic 
approximation, which allows an analytical calculation of the density and temperature profiles |11| . Recently Aumaitre 
et al. (201 investigated, by means of Molecular Dynamics (MD), the fluctuations of the total energy of the system. 
In particular they have studied the behavior of the first two moments of the energy probability distribution function 
(pdf) when the system size is changed, at constant averaged density. Because of the inhomogeneities, the mean kinetic 
energy is no more proportional to the number of particles, and thus it is not an extensive quantity; analogously the 
mean kinetic temperature, defined as the average kinetic energy per particle, is no more intensive. This has led |20) 
to the proposal of a definition of an effective temperature Tg and an effective number of particles Nf , such that Te 
is intensive and the total mean energy (E) = NfTE is extensive in this effective number of particles. In the following 
we will show how an approximate calculation which neglects correlations and the small non-Gaussianities, using the 
hydrodynamic prediction of for the temperature profile, allows in fact to explain the phenomenology observed in 
|20| . Within this description it is possible to obtain an explicit expression of the effective temperature and effective 
number of particles as a function of the system parameters (i.e. number of particles, restitution coefficient, and tem- 
perature of the vibrating wall). At large enough system sizes, the effective temperature is shown to be independent of 
the system size, while the effective number of particles is proportional to the surface of the heating boundary. 



2.1 Energy Probability Distribution Function 

We consider a granular gas of N smooth inelastic hard spheres of diameter a between two parallel walls. The distance 
between the two walls is denoted by H, oriented along the x axis. The linear extension of the wall is denoted by L, 
and periodic boundary conditions are applied along the directions parallel to the wall. One of the walls (in x = 0) is 
vibrated, in order to compensate for the energy lost through inelastic collisions. As usual in the IHS model, the total 
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momentum is conserved in collisions, and only the normal component of the velocity is affected. Thus, the collision 
law for a pair of particles (1,2) reads: 



vj =vi- !(l + a)(vi2 -a)a 
V2 = v 2 + i(l + a)(vi2 ■ cr)cr 



(1) 



where <x is a unitary vector along the center of the colliding particles at contact. In the dilute limit, such a system is 
well described by the Boltzmann equation, involving the one-particle distribution function /(r, v,i): 



ft/tr.vi.tj+v- V/(r,vi,t) = J[f\f] . 
Here J[f\f] is the collision integral, which takes into account the inelasticity of the particles: 

7(r,vr,t)/(r,v 2 *V) 



d-l 



dV5 



d«r(vi2 • <x) 



V12 ■CT>0 



/(r,vi,f)/(r, v 2 ,t) 



(2) 



(3) 



where a collision transforms the couple (v**, v**) into (v*, Vj). The hydrodynamic fields (local density n, local average 
velocity field u, local temperature T) are defined as the velocity moments: 



n(r,t) = J dv/(r,v,t), 
n(r, t)u(r, t) = J dv v/(r, v, t) , 
\n{r,t)T(T,t) = J dv^(v-u) 2 /(r,v,i). 



(4a) 
(4b) 
(4c) 



where d is the space dimension. The hydrodynamic balance equations for the afore-defined fields are derived taking 
the velocity moments in Then, it is possible to show that a steady state solution without macroscopic velocity 
flow exists: this solution is invariant for translations along all directions perpendicular to x, so that the hydrodynamic 
fields only depend upon x. Moreover this solution is stable if the size (in the directions perpendicular to x) of the 
system is smaller than a critical size which depends upon a. In the rest of this section we consider such a "horizontally 
homogeneous" solution, whose temperature profile reads |11| : 



T{1) = T 



cosh ( y/a(a)(i m - i) 



cosh ( yj a(a)£ n 



(5) 



where a(a) is a function of the restitution coefficient (its complete expression is given in ref. jllp. The variable £ is 
proportional to the integrated density of the system over the x axis. Its definition is given by the following relation 
involving the local mean- free-path \{x): 



dx 



X(x) 



r-r d-l 

V2tt— 

r[(d + i)/2]' 



-1 -1 



d-l 



n(x) 



with I = at x = 0. On the other hand, £ m is the maximal lvalue, reached at x = H: 



Ca 



d-l 



dx n{x) = Ca^K 







(6) 



(7) 



where N x is the number of particles per unit section perpendicular to the x axis (e.g. N x = N/L for d = 2), and C 
is a constant depending on the dimension (e.g. C = 2y2 for d = 2). It thus turns out that the crucial quantity £ m 
which determines the hydrodynamic profiles is proportional to the number of layers of particles at rest. The boundary 
conditions used to get expression (|SJl are: 



T(l - 0) = T , 



dT 
~d£ 



e=e„ 



dT 

dx 



, 



(8) 



-H 



where we have assumed that the vibrating wall acts as a thermostat that fixes to To the temperature at x = 0. A 
detailed observation would show that this is not a very realistic assumption. Moreover, the region in the vicinity of 
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x ~ is the one in which the hydrodynamic description may break down. Nevertheless, Tq may be seen as an effective 
parameter such that coincides with the observed profile in the region of the system where hydrodynamics holds. 
In the following we will suppose the velocity distribution to be a Maxwellian at each fixed distance from the wall 
(non-Gaussian features may be observed, but are not relevant for the present calculation), with a local temperature 
(variance) given by JSJl: 

The distribution for the energy of one particle (e = v 2 /2) is hence: 



p(e,£) = f^ i (e) , (10) 



where f a ,v(x) is the gamma distribution |21| : 



Mx) = j^r^- 1 e- ax . (11) 



The characteristic function of the gamma distribution is 

r-+oo 



UAk) = I dxe ikx f a , v (x) = (!--) • (12) 



Our interest goes to the macroscopic fluctuations integrated over all the system. Thus, the macroscopic variable of 
interest is the granular temperature T g , defined here as the average of the local temperature over the x profile: 



T g = I „(r)T(r) dr = f- / " 7 



1 . 1 



To obtain an expression of the energy pdf over all the system, it is useful to divide the box in £ m /A£ boxes of equal 
height (in the I scale) Al, The use of this length scale allows to have a fixed number of particle Ne in each box of size 
L x A£, since dl oc n{x)dx. Moreover, in each box i the temperature can be taken as uniform, denoted Ti = T(i Al). 
In the j — th box, we can calculate the energy ej which is the sum of the energies of the Ne particles in the box. The 
pdf of this energy, qj(x) = prob(e_, = x), can be straightforwardly computed when the velocities of the particles arc 
supposed to be uncorrelated: 

qj{x) = f±. Nl &(x) , (14) 

Thus, the characteristic function for the total kinetic energy E = ^2 e j can be obtained as the product of the charac- 
teristic function qj(k) of the pdf of each tf 

P(k)= TJ q j {k)= n (l-iATi) - ^ ■ (15) 

Since the number of particle per box Ne is a known fraction of the total number of particles (JV^ = NAl/l m ), one can 
rewrite the expression ()15f) as a Riemann sum. In the limit Al — ► this yields the total kinetic energy characteristic 
function: 

P(k) = exp [-^f^ (1 - ikT(£)) uj . (16) 

Note that this result is valid for any temperature profile T(l) and hence it can be applied also to other situations with 
different boundary conditions or different hydrodynamic equations. 



2.2 Local energy fluctuations 

While we have focused on the pdf of the global energy in the previous subsection, it is also interesting to consider an 
intermediate scale and to compute the energy pdf for a " slice" of the system at a given distance x from the vibrating 
wall, and of small height Ax (such that the local density and temperature can be considered as uniform). We still 
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1. Temperature profile in the I scale for a system of 
1000 particles in a box of height H — 200 a and den- 



sity no — 0.05 with a = 0.99. The temperature is scaled 
with the temperature near the vibrating wall To. The in- 
set show the temperature profile of the x and y component 
in the x scale, showing that the temperature is effectively 
isotropic apart from a small region near the vibrating wall. 
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Fig. 2. Probability density function of the y component 
of the total energy E y = Vy. (dots). The parameters of 
the system are the same as in Fig. The solid line is the 
numerical inverse Fourier transform of Eq. (1161 . while the 
dashed line is a \ 2 distribution with same mean and same 
variance of the numerical data, and the dotted line a \ 2 
distribution with parameters corresponding to the granular 
temperature and the true number of degrees of freedom. 
The energy is in units of the temperature of the wall To. 
The inset shows the kurtosis excess Ki = {vf)/(vi ) 2 — 3 for 
the x and y component of the velocity pdf. 



consider that the velocities are Gaussian and uncorrelated. In this case the energy e is a sum of N(x) identically 
distributed random variables, but their number N(x) fluctuates in time. Hence the pdf of the local energy is: 

oo 

Pl(e,x)= Pd(N(x))f^ N{x)i (e)+P d (0)6(e) , (17) 

N(x) = l 

where Pd is the distribution of the number of particles at a distance x from the heated wall. The second term in the 
above expression is due to the fact that when there are no particles the energy pdf is a Dirac distribution. If the 
particles are uncorrelated the distribution Pd should be well approximated by a Poisson distribution with average 
(N(x)) = n(x)Ax: 



P d (N(x) = i) = e -n(«)^ mp _ (18) 



(n(x)Ax) 



i 



In this framework the explicit expression for the probability of the local energy is given, in two dimensions, by: 



Pi(e,x) = e ^y__7 1 ^___J +2e 5(e), (19) 

where /„ (x) is the first kind modified Bcssel function [221 • For higher dimensions the sum in Eq. (|17fl with Poissonian 
density fluctuations can be expressed as a combination of generalized hypergeometric functions. The cumulants of 
those distributions are, in general dimension d: 

(e p )c(x) = "^f^W + 2) . . • (d + 2p - 2) , (20) 
where the notation (X) c denotes the cumulant of the variable X. 



2.3 Comparison with simulations 



To obtain the analytic form of the energy pdf one should calculate the inverse Fourier Transform of equation (|16fl 
using (0) as temperature profile. This docs not seem possible analytically, but one can resort to numerical procedures. 
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Aumaitrc et al. |2HI2()| have shown by Molecular dynamic simulations that the energy pdf is very well fitted by a 
\ 2 law, but with a number of degrees of freedom different from N d, and a temperature different from the granular 
temperature T g . They proposed to define an effective number of degrees of freedom Nf and an effective temperature 
Te, such that energy pdf should be of the form 11(E) — f 1 N f (E), with: 

w '= 2 ffc- T °=%t- (21) 

In Fig.^a direct comparison between the measured temperature profile in the f -scale and the theoretical prediction of 
Eq. (jSJ) shows the agreement between the Molecular Dynamics (MD) simulation and the hydrodynamic prediction. In 
Fig. |21 the Inverse Fourier Transform of (|16fl is compared with the pdf of the y component of the total energy obtained 
from MD simulations, showing a good agreement. We restrict our analysis to the y component because the velocity pdf 
in this direction is closer to the Gaussian distribution (as quantified in the inset by the excess kurtosis, see also |llp. 
and is therefore in better agreement with our hypotheses. In the same figure we have also plotted the function 11(E) 
previously defined, which has a similar shape and yields a better fit of the numerical data. We emphasize however 
that our prediction needs only one fitting parameter, To, which essentially fixes the average of the distribution, while 
in order to have a closed form for 11(E) one needs to know both the mean and the variance of the distribution. This 
underlines the importance of the inhomogencous hydrodynamic profile in determining the fluctuations of the global 
energy; on the other hand, the correlations, which have been neglected in our approximate approach, seem to play a 
lesser role. 

Our approach also allows to investigate the behavior of the macroscopic quantities Nf and Te when the system size 
is changed. Indeed, an expression of the cumulants of the total kinetic energy can be obtained from the characteristic 
function (|16[1 : 

( EP )c = ^~ r^(i)d£ , (22) 



2/ 



so that 

Nd 



(f-T{l)di) t-T\£)dl , 

Nf = , T E = ^ — , 23 

£ m f™ T >{l)dt J^T(£)d£ 

which shows that Te is not rigorously intensive but depends on the number of layers of particles through £ m . At large 
enough £ m however, the integral in equation l|22(l becomes size independent: 

f m T v(t)de~—^= (24) 
Jo 2p^/a(a) 

so that the effective temperature Te defined above becomes a constant proportional to the temperature of the wall, 
while the parameter Nf still depends on the system size: 

A^^f, (25, 
y/a(a) t-m £ 

We now investigate how the numerical results of |2()| fit with the previous framework, for the various possibilities 
of system size variations at fixed particle size a and constant density p = N J (HL) (we consider for simplicity a 
two-dimensional system, as in |2()|): 

— if only L is increased, at constant H, N/L is constant so that £ m and thus Te do not change 

— if only H is increased, at constant T, £ m oc N/L increases so that Te is not constant. The variation of Te with H 
is however slow and moreover, as soon as N is large enough, Te reaches its asymptotic value To/2. 

— if both H and L are increased at the same rate, £ m increases slower so that once again Te is not rigorously constant 
but varies slowly towards To/2. 

The behavior of Nf also depends on the maximum of the integrated density l m . At constant a and density 
p = N/(HL), Eq. (|25() shows that (i) for a square cell, since l m oc \[N , one obtains Nf oc V^V, close to the 
A^ 4 behavior numerically observed in |2D], and (ii) at constant L, if only the height H of the cell is increased, 
£ m is proportional to N, and Nf thus becomes constant. All these features are in agreement with the numerical 
observations in |20| . The above results show that the issues of intensivity and extensivity can be understood through an 
approximate approach which neglects correlations but focuses on the inhomogeneities as described by the hydrodynamic 
profiles |11|. The crucial quantity turns out to be the number of layers of particles at rest. If this number is large 
enough, the effective temperature Te becomes intensive. Moreover, this approach is able to quantitatively describe the 
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Fig. 3. Plot of the effective density rif(x) divided by 
the true density n(x) in a square box at constant density 
pa 2 = 0.04 for a = 0.99 and several system sizes. The 
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Fig. 4. Plot of the local energy pdf for a system of N = 
100 particles in a square cell at density pa 2 = 0.04 and 
a — 0.95. The energy has been sampled in a strip of width 
Ax = 1 .66 o centered in x = H/2. The solid line is the 
analytical estimation from Eq. (IS) . 



behavior of the fluctuations of the total kinetic energy of a vibrated granular gas. In some cases the energy pdf can 
be approximated with a gamma distribution, which is the standard distribution for the energy pdf in the canonical 
equilibrium. Nevertheless there arc strong deviations from the equilibrium theory of fluctuations, since in this case the 
two parameters of the gamma distribution (i.e. the temperature and the number of degrees of freedom) are not the 
granular temperature nor the total number of degrees of freedom. Another important remark is that correlations, and 
in particular contributions from the two points distribution function, do not play a primary role in explaining those 
deviations from the equilibrium theory of fluctuations. In order to characterize deviations which do not arise from 
inhomogeneities, it can be useful to measure energy fluctuations at a given height x from the vibrating wall. From 
the expression of the cumulants of the local energy pdf (|20J) it is possible to define an "effective density" (of particles) 

rif(x) = rf ^1" 2 ' ) ct; ■ If the uncorrelated and Gaussian hypotheses are satisfied, then rif(x) = n(x), and deviations from 
this identity should be due only to correlations and non-Gaussianities. In Figs.|3and0]results from MD simulations are 
shown at various values of the inelasticity, density and number of particles. The agreement between Eq. (|19|) and the 
numerical data is very good, while the deviations from rif(x) — n(x) shown in Fig. El are evident but small compared 
to the deviations between Nf and N. 

Another way to discriminate the respective role of heterogeneities and of correlations for the global system is to 
study granular gases maintained in an homogeneous state by a suitable energy injection, as proposed below. 



3 The homogeneously driven case 

In this section we consider a granular gas kept in a stationary state by an external homogeneous thermostat. Two 
different kinds of thermostat will be used. We will first consider the so called "Stochastic thermostat" , which is a white 
noise acting independently on each particle 24 25 26 27 28 29 30 . In a second part we will show some numerical results 
concerning the "Gaussian thermostat" , which consists of a negative friction force acting on each particle |28l3()j . 



3.1 Stochastic thermostat 

We consider here a gas of N inelastic smooth hard spheres driven by an external random noise. The equation of motion 
of each particle i is: 

^T = -+l(t), (26) 
dr m 

where is the force due to inelastic collisions [the latter being ruled by Eqs. (TJ], and £ 4 is the random acceleration 
due to the stochastic force, which is assumed to be an uncorrelated Gaussian white noise: 

= ^ijS a0 S(t - , (27) 
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line shows a gamma distribution with same mean and same ing uncorrelated velocities 13711 is plotted in dashed line, 
variance. 



where i and j denote the labels of the particles, while a and (3 refer to the Euclidean component of the noise. 
Our interest goes to the stationary state of such system, which is homogeneous, and furthermore does not develop 
spatial instabilities. A good description of such a system can therefore be obtained through a homogeneous Boltzmann 
equation, which will contain a Fokkcr-Planck diffusion term, taking into account the stochastic force. This equation 
reads [25| : 

£o (J_\ 
2 \dvj 

where J[/|/] is the collision integral (cf. Eq. O ■ The stationary solution of this equation is well approximated by the 
product of a Gaussian with a Soninc Polynomial 



&/(vi,t) = J[f\f] 



f(yx,t) 



(28) 



/(c) 



,d/2 



(1 + a 2 S 2 (c 2 )) , 



(29) 



where c = v/vq is a dimensionless velocity, with vq = v2T, S 2 (x) is the second Sonine polynomial 

1 



S 2 (x) 



\{d + 2)x+\d{d + 2) 
2 o 



and a 2 is a coefficient proportional to the kurtosis of the function /(c) 



d(d + 2) 



-]d(d+2) 



An approximate expression of the coefficient a 2 for an arbitrary restitution coefficient is |25j : 

16(1 - a)(l - 2a 2 ) 



a 2 (a) 



73 + 56d - 24ad - 105a + 30(1 - a)a 2 



(30) 



(31) 



(32) 



The energy dissipated in the inelastic collisions is compensated by the energy injected by the thermostat so that 
the system reaches a stationary state, in which the temperature fluctuates around its mean value 



T„ 



(1 — a 2 )f2dna 



d-l 



2/3 



(i + 0K>) , 



(33) 



where f2d is the d-dimcnsional solid angle, m the mass of the particle, and n the density of the system. Here we are 
interested in the fluctuations of the total energy measured by the quantity 



{E*{t)) - {E{t)f 
(E(t)) 2 



(34) 
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Note that a 2 E = 2N/Nf, and that T g is intensive while the total energy is extensive. Brey et al. have computed, by 
means of kinetic equations, an analytical expression for u 2 E in the homogeneous cooling state, which is equivalent to the 
so-called Gaussian deterministic thermostat |14| . This quantity has also been computed in a one dimensional granular 
gas with a similar thermostat jjH]. One of the main differences of the stochastic thermostat with a deterministic one 
is found in the elastic limit. On the one hand, for the cooling state, when the restitution coefficient tends to 1, the 
conservation of energy imposes that the energy pdf is a Dirac delta function, and the quantity <te goes to 0. On the 
other hand, with the stochastic thermostat, if the elastic limit is taken keeping the temperature constant, the strength 
of the white noise will tend to zero, but it will still play a role in the velocity correlation function (we note that it is 
not necessary to keep the temperature constant while performing the limit, since this quantity only sets an irrelevant 
time scale). 

An efficient way to have numerical solutions of Boltzmann-likc equations is the so called Direct Simulation Monte 
Carlo (DSMC). This algorithm simulates a Markov chain whose associated master equation is expected to converge to 
the Boltzmann equation. We performed this kind of simulation to measure the energy pdf. A plot of this quantity is 
shown in Fig. [5] We observe again that it is close to a ^-distribution with same mean and same variance. Nevertheless 
the number of degrees of freedom of this ^-distribution is lower than the true number of degrees of freedom (i.e. 
(N — 1 ) x d). This effect may arise from two separated causes: the non-gaussianity of the velocity pdf, and the presence 
of correlations between the velocities. This feature also suggests that a calculation of the energy pdf with the hypothesis 
of uncorrelated velocities (but non-Gaussian) could explain at least a part of this non-trivial effect. The calculation of 
the pdf of the sum of the square of n variables distributed following 1|29[) is straightforward. The characteristic function 
of the energy pdf is: 

~ /1N f / d(d + 2) ( 1 2 \\ W 

where TV is the number of particles of the system. This results yields 

{E) = %NT, (E 2 )-(E) 2 = ^NT 2 (l + ^a 2 ) , (36) 



2 2 

so that the explicit expression for the energy fluctuations reads: 

, 2 / d+2 \ , . 

> = d{ 1 + — a >) • (37) 

As expected, the temperature does not appear in the dimcnsionless a E above, which depends on the only available 
dimensionless parameters a and d. This result is compared in Fig. with the result of DSMC simulations, carried 
out for several values of the restitution coefficient a and for two different values for the number of particles N. The 
disagreement between the uncorrelated calculation and the simulations is a clear sign of the correlations induced by 
the inelasticity of the system (the tail of the velocity pdf is also non-Gaussian [23 , but this feature has negligible 
consequences for the quantities of interest here). One can note that the fluctuations increase when the restitution 
coefficient decreases. One can also see that there is a value of the restitution coefficient a around that is when 

the approximate expression of a 2 vanishes, for which a E is exactly f = 2/d, as for a gas in the canonical equilibrium. 

We now turn to the dependence of a E on the strength of the white noise r. It is useful, for this purpose, to 
introduce a rcscaled, dimcnsionless energy 

E = E ^ E) . (38) 

We have plotted in Fig. [7| this rescaled energy pdf for a system of N = 100 particles with a restitution coefficient 
a = 0.5 and for several values of the strength of the white noise r. As expected, all the pdfs collapse into a unique 
distribution: the role of the noise's strength is only to set the temperature (or mean kinetic energy) scale. Besides, 
relative energy fluctuations depend only on a and N. Moreover, since a 2 E does not depend on the number of particles 
N (for N large enough), J (E 2 ) — (E) 2 grows as y/N, so that the central limit theorem applies, and hence P(E) is 
a Gaussian in the thermodynamic (N — ► oo) limit. Figure |H1 indeed shows how the rescaled energy pdf for N = 1000 
particles is very close to a Gaussian, even if the x 2 -distribution still gives a better fit. 



3.2 Gaussian thermostat 

In this section, we will focus on the fluctuations in a granular gas driven by a Gaussian thermostat , which consist in 
adding a viscous force, with a negative friction coefficient, in the equation of motion of each particle. Thus the velocity 

1 This is the traditional nomenclature to denote such a thermostat for granular gases. Note, however, that the thermostating 
force acting on the particles is not obtained in such a way to satisfy Gauss' principle of least constraint. The total kinetic energy 
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Fig. 7. Plot of the pdf of E for a restitution coefficient 
a — 0.5, for N = 100 particles, and for several values of the 
noise's strength The simulations have been carried on 
with both Monte-Carlo (DSMC) and Molecular Dynamics 
(MD) algorithms, with na 2 = 0.04. 
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8. Plot of the pdf of E for a restitution coefficient 
a = 0.5, for TV = 1000 particles, and for several values 
of the noise's strength £q. The solid line is a rescaled \ 2 ' 
distribution with Nf = 2N/a% ~ 1759 degrees of freedom, 
and it is very close to the Gaussian distribution, plotted 
here in dashed line. 



of the i-th particle will verify the following equation: 

^ = 5+7v 4 . (39) 
at m 

In the above equation F c is the force due to collisions, and 7 is a positive constant. The Boltzmann equation corre- 
sponding to Eq. is: 

&/(vi,t) = - 7 A • [ Vl /( Vl ,i)] . (40) 

In this case also, the relevant solution is well approximated by a Gaussian multiplied with a Sonine polynomial [cf. 
Eq. (|29|l ]. The expression of the coefficient 02 has been widely studied [25128133] . and the one which is in better 
agreement with numerical simulations is the following |28j : 

a2(a) = 16(1-^(1^) 

2{ ) 25 + 24d-a(57-8d)-2(l-a)a 2 ' 1 ' 

The parameter 7 determines the strength of the homogeneous force acting on the particles. Its primary role is to set 
the temperature scale of the stationary state. Projecting Eq. I|40|l on the second velocity moment is it possible to get 
an approximate expression of the stationary temperature: 

r -( a-^- ) 2(1 + 0M) - (42) 

This result has been obtained assuming Gaussian velocity pdf, since the corrections coming from the coefficient 02 are 
small. 

Equation l|40|l is formally equivalent to the Boltzmann equation for the homogeneous cooling state (HCS) of a 
granular gas, for the particular solution in which all the time dependence of the velocity pdf appears trough the time 
dependent temperature |34| . This mapping of the cooling state onto a steady state can be obtained in an equivalent 
fashion applying a logarithmic time-scale transformation to the time-dependent Boltzmann equation of the HCS |35I14| . 
Furthermore, the above scaling solution can be extended to the n-particle distribution function of the system. The 
resulting hierarchy equations of the HCS are then exactly the same one would obtain for a homogeneous granular 
gas thermostated by a Gaussian thermostat. Thus these two models not only have the same velocity pdf, but also 
the same n-particles distribution, and hence they will have the same relative fluctuations and correlation functions. 
Henceforth we will no more make any distinction between the HCS and the Gaussian thermostat. We will therefore 



is not constant, but fluctuates in time. Henceforth we will refer to the Gaussian thermostat without any connection to Gauss' 
principle of least constraint, nor to the thermostats described in |32j . 
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Fig. 9. Plot of a E versus the restitution coefficient a for 
N = 100 (O) and N = 1000 (□) particles, for the Gaussian 
thermostat. The solid line shows the theoretical predictions 
of [Tl], given by Eq. ijl^jl . 
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Fig. 10. pdf of the rescaled energy E when the restitu- 
tion coefficient is a = 0.5. The system is driven with the 
Gaussian thermostat with 100 (O) an d 1000 (□) particles. 
The solid line shows a rescaled gamma distribution with 
zero- mean, a variance equal to 1, and with a number of 
degrees of freedom Nf given by 2N/a%. The dashed line is 
the Gaussian with zero-mean and unit variance. 



compare results obtained for the HCS with numerical simulations of a gas driven in a stationary state with a Gaussian 
thermostat. We also have to note that this kind of homogeneous state is unstable under small density fluctuations, 
and naturally develops cluster and shear instabilities. In the following we will therefore focus only on homogeneous 
states. For the HCS Brey et al. ^2] have obtained an analytical expression for the quantity a 2 E defined in Eq. 
Its expression is: 

n d(d+l) did + 2) ,,,,,,, , . 

4 = V 2 ' + 1 4 ' a 2 (a) + d 2 b(a) , (43) 

where 2 

4 (1- 2a 2 ) (1- 3a + a 2 (a- 1)) - (1 - a) (3 -2a (8 + 5a)) d- 2 (3 + a) d 2 
6(<l) = 2d(7 + 6d-a (15 + 2a (1- a) -2d)) ' (44) 

Note that 6(1) = —(1 + d)/{2d) so that a\ — > as a — ► 1. The results from DSMC simulations and the theoretical 
predictions by Brey et al. are shown in Fig. |5J The simulations have been carried out for two different sizes of the 
system (N = 100 and N = 1000), and for several values of a, and in each case the numerical results are in good 
agreement with the predictions of the relation l|43|l . The energy pdf behaves again as a \ 2 distribution. When the 
system is not too large (N ~ 100) the shape of the energy pdf is very well fitted with a \ 2 distribution with a 
mean value obtained from the temperature expression (|42|) and a variance obtained from Ij43() . Furthermore, when the 
number of particles is increased, the ^-distribution becomes closer and closer to the Gaussian distribution. All these 
features are shown in Fig. ^] Ref . ^1] reported a good agreement between Molecular Dynamics results and Eq. I|43(l . 
The fact that the relevant correlations are already captured by the DSMC scheme is a remarkable feature that will be 
commented further in the conclusion. 



4 Conclusion 

In this paper, we have investigated the fluctuations of the total kinetic energy of granular gases maintained in a non- 
equilibrium stationary state through various kinds of energy injection mechanisms. For heating through a boundary, 
the system remains inhomogencous, and the characteristic function of the energy pdf can be analytically computed 
as a functional of the temperature profile of the system, under the assumption that the particle velocities are uncor- 
rected. Despite this approximation, the use of the hydrodynamic profiles obtained in |llj allows to recover results in 
agreement with the numerical simulations carried out in Ref. |20j . This result seems to hold even if the hydrodynamic 
approximation is not valid in the whole volume of the system, in particular near the vibrating wall. Moreover, we find 
that the effective temperature defined in [501, which is measurable, is intensive only if the number of layers of particles 

2 Note that the following expression is not exactly the same given in |14|. since here a different expression of the coefficient 
ai (which fits better with numerical simulations) is used. 
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is large enough, and can be related to a quantity which plays the role of the "temperature of the vibrating wall" . 
which is microscopically not well defined, and a priori not measurable. An effective number of degrees of freedom 
can be also defined, and its dependence on the true number of degrees of freedom has been obtained. Once again, 
the obtained behavior is in agreement with the numerical simulations described in and leads to the conclusion 
that the non-extensive behavior of the total energy in such a system is essentially due to the density and temperature 
inhomogeneities and can be understood through a hydrodynamic approach. 

In homogeneously driven granular gases the total energy displays a x 2 pdf, with a number of degrees of freedom 
proportional to the number of particles. The total energy is hence extensive, and moreover, in the N — > oo limit the 
X 2 distribution tends to the Gaussian. We measured, by means of DSMC simulations, the rescaled variance a E of 
the energy pdf, which depends only on the restitution coefficient. These non-trivial fluctuations are essentially due to 
the correlations induced by the inelasticity of the particles. We can distinguish between two different contributions to 
these correlations. First, the non-Gaussianity of the velocity pdf, which simply tells that the Euclidean components of 
the velocity of each particle are correlated one to each other. Second, a contribution from the two particles velocity 
pdf, which does not factorize exactly as a product of two one-particle distributions. The analytical computation of the 
rescaled variance o~ E under the assumption of uncorrelated non-Gaussian individual velocity pdfs does not reproduce 
the DSMC results, showing the relevance of this second contribution. It must be pointed out, however, that these 
correlations do not invalidate the Boltzmann equation. As already noted in [13114] , the two points correlation function 
ff2(vi,V2), which is defined by: 

5 2 (vi,v 2 ) = / (2 >(vi,v 2 ) - /(vi)/(v 2 ) , (45) 

where is the two points distribution, is of higher order in the density expansion (roughly speaking O (g(vi, v 2 )) ~ 
® (/( v i)/( v 2)) /N). This is confirmed by the numerical observation, since the energy pdf tends to a Gaussian when 
the number of particles increases. 

When the gas is heated with the Gaussian thermostat, the numerical results are in good agreement with the 
analytical result in ^3]. This means in particular that the DSMC algorithm is able to capture effects arising from 
the two-particle velocity pdf, which is remarkable. Moreover the energy pdf obtained with DSMC coincides with its 
counterpart found in MD simulations in the dilute limit, provided the system remains homogeneous (cf. Fig. UJ. This 
feature is consistent with the fact that the DSMC algorithm exactly simulates the homogeneous pseudo-Liouvillc 
equation governing the A^-particle dynamics of the gas. 

For the gas thermostated by a random noise, numerical simulations show that the quantity a\ behaves in a fashion 
similar to what happens for the Gaussian thermostat. Nevertheless in the elastic limit (a — ► 1) it tends to a non 
trivial value (~ d/^/2), which is still not understood. An analytical calculation of o~ E is in principle possible, exploiting 
the methods described in Nevertheless, in order to carry on with such methods, one needs the expression of the 
eigenfunctions of the linearized Boltzmann operator. While for the Gaussian thermostat these eigenfunctions can be 
computed, in the way as for the elastic hard sphere gas, this task seems more demanding in the case of a stochastically 
driven granular gas. 

The authors thank S. Aumaitre, S. Fauve and J. Farago for stimulating discussions. 
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